function [ F ] = getDeformationGradientFD( x, y, z, mesh, theta )
%GETDEFORMATIONGRADIENTFD Summary of this function goes here
%   Detailed explanation goes here

epsilon = 0.00001;

nTri = size(mesh.f, 1);

F = zeros(3, 3);

% up

up = zeros(3, 1);
p = [x, y+epsilon, z];

for j = 1 : nTri        
        
    q0 = mesh.v(mesh.f(j, 1), :);
    q1 = mesh.v(mesh.f(j, 2), :);
    q2 = mesh.v(mesh.f(j, 3), :);
              
    T = triarea(q0, q1, q2);
        
    v0 = evaluateQuadraturePoint(p, q0);
    v1 = evaluateQuadraturePoint(p, q1);
    v2 = evaluateQuadraturePoint(p, q2);
        
    up(:) = up(:) + ((v0 + v1 + v2) * T / 3.0 * theta(3*j-2:3*j)); 
end

% down

down = zeros(3, 1);
p = [x, y-epsilon, z];

for j = 1 : nTri        
        
    q0 = mesh.v(mesh.f(j, 1), :);
    q1 = mesh.v(mesh.f(j, 2), :);
    q2 = mesh.v(mesh.f(j, 3), :);
              
    T = triarea(q0, q1, q2);
        
    v0 = evaluateQuadraturePoint(p, q0);
    v1 = evaluateQuadraturePoint(p, q1);
    v2 = evaluateQuadraturePoint(p, q2);
        
    down(:) = down(:) + ((v0 + v1 + v2) * T / 3.0 * theta(3*j-2:3*j)); 
end

% left

left = zeros(3, 1);
p = [x-epsilon, y, z];

for j = 1 : nTri        
        
    q0 = mesh.v(mesh.f(j, 1), :);
    q1 = mesh.v(mesh.f(j, 2), :);
    q2 = mesh.v(mesh.f(j, 3), :);
              
    T = triarea(q0, q1, q2);
        
    v0 = evaluateQuadraturePoint(p, q0);
    v1 = evaluateQuadraturePoint(p, q1);
    v2 = evaluateQuadraturePoint(p, q2);
        
    left(:) = left(:) + ((v0 + v1 + v2) * T / 3.0 * theta(3*j-2:3*j)); 
end

% right

right = zeros(3, 1);
p = [x+epsilon, y, z];

for j = 1 : nTri        
        
    q0 = mesh.v(mesh.f(j, 1), :);
    q1 = mesh.v(mesh.f(j, 2), :);
    q2 = mesh.v(mesh.f(j, 3), :);
              
    T = triarea(q0, q1, q2);
        
    v0 = evaluateQuadraturePoint(p, q0);
    v1 = evaluateQuadraturePoint(p, q1);
    v2 = evaluateQuadraturePoint(p, q2);
        
    right(:) = right(:) + ((v0 + v1 + v2) * T / 3.0 * theta(3*j-2:3*j)); 
end

% near

near = zeros(3, 1);
p = [x, y, z-epsilon];

for j = 1 : nTri        
        
    q0 = mesh.v(mesh.f(j, 1), :);
    q1 = mesh.v(mesh.f(j, 2), :);
    q2 = mesh.v(mesh.f(j, 3), :);
              
    T = triarea(q0, q1, q2);
        
    v0 = evaluateQuadraturePoint(p, q0);
    v1 = evaluateQuadraturePoint(p, q1);
    v2 = evaluateQuadraturePoint(p, q2);
        
    near(:) = near(:) + ((v0 + v1 + v2) * T / 3.0 * theta(3*j-2:3*j)); 
end

% far

far = zeros(3, 1);
p = [x, y, z+epsilon];

for j = 1 : nTri        
        
    q0 = mesh.v(mesh.f(j, 1), :);
    q1 = mesh.v(mesh.f(j, 2), :);
    q2 = mesh.v(mesh.f(j, 3), :);
              
    T = triarea(q0, q1, q2);
        
    v0 = evaluateQuadraturePoint(p, q0);
    v1 = evaluateQuadraturePoint(p, q1);
    v2 = evaluateQuadraturePoint(p, q2);
        
    far(:) = far(:) + ((v0 + v1 + v2) * T / 3.0 * theta(3*j-2:3*j)); 
end

px = (right - left) / (2 * epsilon);
py = (up - down) / (2 * epsilon);
pz = (far - near) / (2 * epsilon);

F(1,1) = 2 * px(1);
F(2,2) = 2 * py(2);
F(3,3) = 2 * pz(3);
% F(1,2) = px(2);
% F(1,3) = px(3);
% F(2,3) = py(3);
% F(2,1) = py(1);
% F(3,1) = pz(1);
% F(3,2) = pz(2);
F(1,2) = px(2) + py(1);
F(1,3) = px(3) + pz(1);
F(2,3) = py(3) + pz(2);
F(2,1) = F(1,2);
F(3,1) = F(1,3);
F(3,2) = F(2,3);

end

